library(list)
library(readxl)

soc.all <- read.csv("Data/data_list_experiment_for_R_final.csv", header = TRUE, sep = ",", fileEncoding = "latin1")

soc <- subset(soc.all, dum_2==0 & dum_3==0 & sample_part==1 & Referer!=3)
soc.dir <- subset(soc, dum4_1==1)
fit.list <- ictreg(v16_1 ~ 1, treat="dum_1", J=4, data=soc, method="lm")
summary(fit.list)
fit.sens <- glm(v17_1 ~ 1, data=soc.dir, family=binomial("logit"))
summary(fit.sens)
avg.pred.social.desirability1 <- predict(fit.list, direct.glm=fit.sens, family=binomial("logit"))
ceiling.results <- ictreg(v16_1 ~ 1, treat = "dum_1", 
                          J=4, data = soc, method = "ml", fit.start = "nls",
                          ceiling = TRUE, ceiling.fit = "bayesglm")
summary(ceiling.results)

soc <- subset(soc.all, dum_1==0 & dum_3==0 & sample_part==1 & Referer!=3)
soc.dir <- subset(soc, dum4_2==1)
fit.list <- ictreg(v16_1 ~ 1, treat="dum_2", J=4, data=soc, method="lm")
summary(fit.list)
fit.sens <- glm(v17_1 ~ 1, data=soc.dir, family=binomial("logit"))
summary(fit.sens)
avg.pred.social.desirability2 <- predict(fit.list, direct.glm=fit.sens, family=binomial("logit"))
ceiling.results <- ictreg(v16_1 ~ 1, treat = "dum_2", 
                          J=4, data = soc, method = "ml", fit.start = "nls",
                          ceiling = TRUE, ceiling.fit = "bayesglm")
summary(ceiling.results)

soc <- subset(soc.all, dum_1==0 & dum_2==0 & sample_part==1 & Referer!=3)
soc.dir <- subset(soc, dum4_3==1)
fit.list <- ictreg(v16_1 ~ 1, treat="dum_3", J=4, data=soc, method="lm")
summary(fit.list)
fit.sens <- glm(v17_1 ~ 1, data=soc.dir, family=binomial("logit"))
summary(fit.sens)
avg.pred.social.desirability3 <- predict(fit.list, direct.glm=fit.sens, family=binomial("logit"))
ceiling.results <- ictreg(v16_1 ~ 1, treat = "dum_3", 
                          J=4, data = soc, method = "ml", fit.start = "nls",
                          ceiling = TRUE, ceiling.fit = "bayesglm")
summary(ceiling.results)

pdf(file="Exhibits/Figures_Appendix/figure_C2.pdf")
plot(c(avg.pred.social.desirability1,avg.pred.social.desirability3, avg.pred.social.desirability2), axes.ict=TRUE,
     main = "Figure C2: Comparison between List Experiment and Direct Question Estimates",
     cex.main = 0.9,
     ylab="Estimated Proportions / Difference in Proportions", 
     pch=c(19, 19, 19, 17, 17, 17, 15, 15, 15))
axis(1, at=c(2,5,8), lwd=0, lwd.ticks=0, lab=c("Friends were going", "Wanted to tell friends", "Wanted to post on SM"), cex.axis=0.8)
text(1, 0.2, "List", cex=0.8)
text(2, 0.24, "Direct", cex=0.8)
text(3, 0.24, "Difference", cex=0.8)
text(3, 0.2, "List-Direct", cex=0.8)

dev.off()


